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Abstract 

I develop an improved Hamiltonian for classical, Minkowski Yang-Mills theory, 
which evolves infrared fields with corrections from lattice spacing a beginning at 
0(0"^). I use it to investigate the response of Chern-Simons number to a chemical 
potential, and to compute the maximal Lyapunov exponent. Both quantities have 
small a limits, in both cases within 10% of the limit found using the unimproved 
(Kogut Susskind) Hamiltonian. For the maximal Lyapunov exponent the limits differ 
by about 5%, significant at about 5a, indicating that while a small a limit exists, 
its value is corrupted by lattice artefacts. For the response of Chern-Simons number 
the statistics are not good enough to resolve 5% differences, but it seems possible 
in analogy with the Lyapunov exponent that the final answer depends on the lattice 
regulation. 

1 Introduction 

Several interesting and outstanding questions in high energy physics and cosmology, 
such as the production of quark-gluon plasma in relativistic heavy ion collisions and 
the conditions in the early universe under which baryogenesis may have occurred, 
involve finite temperature field theories where there are interactions between bosons, 
and some of the bosons have very small masses in comparison to the temperature, 
m -C T. The infrared physics of such systems can be strongly coupled; in per- 
turbation theory, the contribution of loops involving very soft {0{g^T), g a generic 
coupling strength characterizing the theory) momenta can be order 1, even when g 
is small enough that the vacuum theory is weakly coupled. This strongly coupled 
infrared physics is often very interesting; in the electroweak theory it may determine 
the strength of the electroweak phase transition, and it sets the rate of Chern-Simons 
number {Ncs) diffusion, which determines the rate of baryon number violation. The 
former problem is thermodynamical, and there are Euclidean tools which can be 
brought to bear on the nonperturbative physics involved[|l], ^. However, the latter 
question is dynamical, and cannot be addressed by Euclidean methods. There is also 
interest in a more general understanding of the infrared dynamics of Yang-Mills and 
Yang-Mills Higgs theory at finite temperature, for instance to understand its impli- 
cations in heavy ion collisions. These questions cannot be addressed by Euclidean 
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methods, which are only appropriate for answering thermodynamical questions. They 
are also very hard to address analytically in any reliable way, and the prospects for 
numerical evaluation of Minkowski path integrals are more or less hopeless. 

Several years ago, Grigoriev and Rubakov made an interesting proposal for a new 
numerical approach to the understanding of the infrared problem in finite temperature 
bosonic systems^. The key observation is that, if the coupling constant g is small, 
then the nonperturbative physics arises because occupation numbers are large in 
the infrared; but this is precisely the regime where a classical field approximation 
becomes reasonable. The physics is nonperturbative but in an essentially classical 
way. Therefore, it is reasonable to hope that, by analyzing the analogous classical 
field theory at finite temperature, one can understand at leading order the physics of 
the quantum system of interest. 

Based on this idea, there has been a revival of interest in classical, thermal Yang- 
Mills theory ^, H, ^ H, ^. The work is based on lattice implementation of clas- 
sical Yang-Mills theory, using the Hamiltonian formulism worked out by Kogut and 
Susskind years ago|lTO|. This work has yielded interesting results; it appears that 
we now know the Ncs diffusion rate in the symmetric electroweak phase to within 
around 5% accuracy p|, and there is numerical evidence for the interesting conjecture 
that the maximal Lyapunov exponent of Yang- Mills theory is twice the damping rate 
of a plasmon at rest0. 

However, there have been some frustrations in the application of this technology; 
for instance, in I find that, to determine the rate of Ncs diffusion in SU(3) gauge 
theory, the lattice spacing must be so fine that the calculation is computationally 
prohibitive. Also, there is a more fundamental challenge to the whole approach. 



raised by Bodecker et. al. |]rT[, who consider the "hard thermal loops," the one 
loop contributions of ultraviolet modes to the dispersion relations and interactions 
of the infrared modes. On the lattice, the strength of the hard thermal loops grows 
linearly with inverse lattice spacing, and as is often the case with linearly divergent 
quantities, the functional form of the hard thermal loops may depend on the form of 
the cutoff; they posess lattice artefacts which make their q/q^ dependence different 
than in the continuum, quantum theory. Although the results for A'';^^ diffusion in 
SU(2) apparently have a fine lattice limit [§], which must be independent of the overall 
amplitude of the hard thermal loops, the value of the fine lattice limit may depend 
on the functional form of the hard thermal loops, and hence on the specifics of the 
lattice cutoff. I brieffy explored this possibility in |^, where I found some evidence 
for a weak dependence on the form of the cutoff; however the case was far from clear, 
and it is an open and interesting question whether the Ncs diffusion rate depends 
on the presence of hard thermal loops and on their functional form, and whether the 
convergence to a fine lattice limit represents the disappearance of nonrenormalizeable 
operators or a suppression of effects which die as a power of ?Ti|,/mno„pgrturbative- 

One or both of these problems might be solved by finding an "improved" Hamil- 
tonian for the classical, Minkowski Yang-Mills theory. By "improved" I mean a 
Hamiltonian such that the interaction of infrared modes would be the same as in 
the continuum, plus finite a effects starting at a^, rather than at a^, as they do in 
the simplest lattice implementation. That is, for infrared fields the lattice action can 
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be expanded in operator dimension, and an "improved" action is one in which, in 
addition to the desired dimension 4 F^^F^^ term, there are no dimension 6 operators, 
but only those starting at dimension 8. It is well known that in evolving scalar par- 
tial differential equations, and in Euclidean lattice gauge theory, such improvement 
greatly reduces the numerical demands by allowing accurate results on much coarser 
lattices[Q. Implementing such an improved Hamiltonian for the systems under con- 
sideration should either lighten the numerical demands, or verify that hard thermal 
loops are indeed important, either because the lattice fineness requirement arises 
from the need to make the Debye mass large, or because the results depend on the 
functional form of the hard thermal loops, and hence on the form of the regulation. 

The basic idea of an improved lattice Hamiltonian is quite simple, and is easiest 
to present in the case of scalar field theory. Consider the evolution of the continuum 
system derived from the Lagrangian 
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(1) 



with a single component, real scalar field. The Hamilton's equations of motion 
(defining vr as the conjugate momentum of cf)) are 

^ = v20 + m20 + A03. (3) 
at 

To implement these equations on a lattice, we discretize the fields (f) and vr so they 
each take on a single value at each lattice point. The only difficulty lies in choosing 
a lattice definition of V^. If we make the replacement 

d^xiVct^f = a'Y.Y. ^^^^ + (4) 

in the Lagrangian, then the appropriate substitution in the equation of motion is 

(X 

i 

which is the standard lattice Laplacian. 

Now consider the quality of this implementation when evolving some smooth field 
configuration. If the characteristic wavelength of the field configuation of interest is 
on order the lattice length, it is inevitable that the lattice implementation will poorly 
reproduce the continuum dynamics; so it is only interesting to test the quality of the 
lattice implementation for fields which are smooth on the lattice scale, in which case 
we can study the quality of the lattice Laplacian by expanding in a Taylor series 
about the point x; it is straightforward to show that 



cpjx - ai) + 0(x + ai) - 20(x) ^ ^ 
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where I write in terms of a positive 3 dimensional measure, as I will for the remainder 
of the paper. All but the first of the derivative terms are undesirable. The second 
term will produce lattice artefacts in the evolution of a mode of wave vector k which 
will be on order a^fc^; if a <^ which is necessary for the lattice implementation to 
at all resemble the desired continuum system, then the quality of the implementation 
would be improved if we could write a new lattice Laplacian which did not produce 
an O(a^) term. 

This is quite easy; just include contributions from next nearest neighbors as well 
as nearest neighbors, and choose the coefficients so the coefficients of and are 
as desired. The correct choice is 
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2 _ ^ ji<Pix + ai) + (f){x - ai)) - Y^{(f){x + 2ai) + (j){x - 2ai)) - ^(f){x) 
i " 

= EV^0 + Oa2EVf0-^EV^0 + .... (8) 

which is free of O(a^) contamination. We can get this form if we replace the lattice 
Lagrangian term —{(j){x) — (f){x — a(z)))^/2 with 

(0(a;) - (p{x - ai)) + + ~ ~ + ^2^^^ ~ ^"^0 ' ^^"^ 

which means we can still expect to have a conserved energy, which must however be 
defined in terms of the above gradient term. 

This system is more complicated to evolve, perhaps doubling the number of com- 
putations in an update. What we have gained, though, is that the lattice size neces- 
sary to model the desired system with a given accuracy is greatly reduced; the number 
of lattice points needed is order the square root of the number previously required. 
This greatly reduces the numerical demands of the evolution. 

The implementation of an improved lattice action is not quite as simple in non- 
abelian gauge theory, but the techniques involved have been worked out in the Eu- 
clidean case some time ago[|13|. The purpose of this paper is to develop and implement 
the correct improved Minkowski Hamiltonian, and to use it to study the cutoff de- 
pendence of the system. 

The paper is organized as follows. Sec. |^ will review the Kogut-Susskind imple- 
mentation of Minkowski, classical lattice Yang-Mills theory, the lattice definition of 
dNcs/dt, the addition of a chemical potential for Ncs, and the thermalization of the 
lattice system. Everything in the section has appeared in previous literature; I include 
it for didactic reasons, to review necessary tools and ideas for the remainder of the 
paper. In Sec. ^I will construct a lattice Hamiltonian in continuum time which is free 
of dimension 6 lattice artefacts and find a discrete time algorithm for its evolution. 
The kinetic (electric field) part of the Hamiltonian is not diagonal in the obvious 
basis for electric fields, so the modifications to the definition of dNcs/dt, chemical 
potential, and the thermalization algorithm are nontrivial. With the Hamiltonian, 
chemical potential, and thermalization fully implemented, I analyze the motion of 
Ncs under a chemical potential in SU(2) gauge theory, presenting numerical results 
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in Sec. |. I also compute the maximal Lyapunov exponent for the system. Finally, 
Sec. 1^ concludes. 



2 Kogut-Susskind Hamiltonian and Ncs 

In this section I will review the construction of a lattice, Minkowski model for classical 
Yang- Mills theory, developed by Kogut and Susskind [|10]- I present a leapfrog algo- 
rithm for its update and a definition for ANcs/^t due to Ambjorn et. al. 0, |^, and 
a way to apply a chemical potential for Ncs developed in 0. Krasnitz has developed 
a thermalization algorithm for the system, based on a set of Langevin equations [0; 
but I will present a much simpler thermalization algorithm developed in 0, because 
it is this algorithm I will extend to the improved model in Sec. ^. 



2.1 Degrees of Freedom and Hamiltonian 

To implement Yang-Mills theory on a lattice, Kogut and Susskind followed Wilson's 
ideas JI^ , taking very literally the nature of the guage field as an affine connection for 
SU(2) objects in space. If there were a scalar field taking values in the fundamental 
representation of SU(2) at each point on a lattice, the gauge field should provide a rule 
for parallel transporting the field from lattice site to lattice site. In the continuum 
the rule is that the value of (j){x + ai), parallel transported to the point x, is 

0^(x + ai)=Pexp / _ —Ai{y)dyj\ (j){x + ai) , (10) 

\J x+ai Z J 

with the path ordering placing points near x to the left. (Here and throughout I will 
use nonrelativistic notation and a positive metric on the space indicies.) The lattice 
analog is to have a matrix \J gSU(2) which lives on the link between the lattice sites, 
so the parallel transported object is just \Ji{x)(^{x^ai). The left index of \J lies at the 
basepoint of the link and the right index lies at its endpoint, so to perform a gauge 
transform in which 0(a;) is rotated by G gSU(2) to the link Vi[x) is rotated to 
GUi{x) and the link Ui{x — ai) becomes Ui{x — ai)G^] so transporting G(j){x) back to 
X — ai gives Ui{x — ai)G'^G(f){x), which is unchanged, as it should be, by the guage 
transform at x. To forward transport 4>{x) to x + ai we multiply by Uj{x). Note that 
my convention is to label the link by its basepoint, so Ui{x) runs from x to x + ai; 
also from now on I will drop a, measuring everything in lattice units, exept in those 
cases where it is useful to include it to count dimensions. I will also drop hats where 
the meaning is obvious, ie x + i means x + ai. 

For the time being I will consider a spatial lattice but with continuous time. The 
connection in the time direction then remains a Lie algebra element living at each 
lattice point, with parallel transport of a field (j){x, t) to time t + At given by 

/ rt+At \ 

(j)t+At{x,t) = Pexp / zt'^A^^ (j){x,t) (11) 
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where the path ordering places earUer times to the right. Aq transforms under a gauge 
change as Aq G{x, t)AQ{x, t)G'^{x, t) + doG^x, t). 

The field strength is a curvature tensor and it refiects the failure of an object, 
on transport around a closed curve, to return to its original value. The smallest 
nontrivial closed curve is the square plaquette, so in the continuum case 

corresponds to the product of the links around a plaquette, 

Uj{x)Ui{x + j)U]{x + i)Ul{x) = Uu{x) (12) 

not being the identity. A reasonable definition of the magnetic field strength is the 
Lie algebra element 

^TV-iT"[/D(x), (13) 

which, although it transforms under gauge change as an adjoint object at the point 
X, is really associated with the plaquette lying in the positive i and j directions from 
X. The continuum Lagrangian depends on the F^j as / F^jFij/A, which is reproduced 
at leading order of a weak field expansion by 

-Lu^J^l-^TrUu (14) 

where the sum is over all plaquettes in the lattice. (The coefficient does not agree 
with the continuum one, but I will account for this when I thcrmalizc the system.) 

Next consider the electric field part of the Lagrangian. It is convenient to write 
the time derivative of a link variable as 

= ir'^mxnix) (15) 

and to expand the Lie algebra element C/f (x) as 

il^ix) = Ao{x) - Ui{x)Ao{x + i)Uj{x) + E^{x) . (16) 

The field transforms under gauge transformations as an adjoint object and can be 
recognized in the weak field limit as the electric field; alternately we can notice that 
it is the small At limit of a plaquette in the i and time direction, traced against the 
Lie algebra and divided by At. The electric part of the Lagrangian should be 

l.=2:™m. (17) 

x,i 

Now we can derive Lagrange equations of motion. Note first that the time deriva- 
tive of Aq never appears, so Aq enters the Lagrangian like a Lagrange multiplier. 
Varying with respect to Aq gives 

= Efe(^)-4(^--^)^Tr-^r"t/j(a;-j)^/[/,(a;-j)) , (18) 



dA${x) 
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where the second term is just the E field on the hnk leading into x, parallel transported 
so its indicies reside at the point x. This equation is a constraint, called the Gauss 
constraint because it corresponds to the continuum condition that DjEj = 0, which 
is Gauss' Law. If we initialize the evolution by choosing values for the U and E, then 
the Gauss constraint is a restriction on the allowed initial conditions. 
The equations of motion derived by varying with respect to a link are 

-^ = -DFnx), DF^x) ^ - ^r^Un , (19) 

where 4nj(x) refers to the 4 plaquettes which contain the link leaving x in the i 
direction, with orientation chosen so that Ui{x) and not Ujlx) appears, and beginning 
and ending at the point x. I use the notation DF°-{x) because it is compact and 
remeniscient of the equivalent continuum quantity, {DjFji)"-{x). Note that Eq. (|T^ 
automatically preserves the Gauss constraint if the initial conditions satisfy it, because 
the time derivative of Eq. (|TB]) will, according to Eq. (0), consist of the sum of 
traces over plaquettes, and every plaquette which contributes to one E contributes 
to a neighbor with opposite sign. This exact satisfaction is necessary because the 
evolution would not make any sense if it excited the modes constrained to be zero; in 
effect we keep track of more E fields than are actually dynamical degrees of freedom, 
and we must make sure that only the dynamical degrees of freedom actually get 
excited. 

Note that Eq. (|T^ does not actually uniquely specify the evolution of the system. 
E!^{x) does not give us and the equations say nothing about how to evolve the 

value of Aq. This is because the evolution should not in fact be uniquely specified, yet; 
we have the freedom to make arbitrary time dependent gauge transformations. To 
write a numerical update algorithm we have to decide how to handle this freedom; the 
simplest and easiest choice is to pick the gauge such that Aq is everywhere zero. The 
system can then be considered as a constrained Hamiltonian system with coordinates 
U and conjugate momenta E. The Hamiltonian is just Le — Lu = He + Hu. 



2.2 Numerical evolution of system 

One way to evolve this system is to use a leapfrog algorithm in which the U fields are 
defined at integer multiples of At, and the E are taken to be Lie algebra elements 
living at half time steps; the rule is 

Ef{x, t + At/ 2) = E^{x, t - At/2) - AtDF^{x, t) (20) 
Ui{x,t + At) = exp{iT^E^{x)At)Ui{x) . (21) 

This choice still identically preserves the Gauss constraint whether or not the 
parallel transport involved in the definition of the Gauss constraint is performed with 
the past or future U. The reason is that E commutes with its own exponent, so the 
parallel transport is the same for the updated U as for the U before update. This 
algorithm forms the basis of the work in 0, |, ^. 
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An alternate evolution algorithm developed in 0], which is closer to what I will 
need for the improved Hamiltonian system, is to begin with a discrete time lattice, 
so the Aq fields are also connection matricies and the E fields are plaquette matricies 
for plaquettes with one space and one time step. The action is 



1 - iTrf/^ 



(22) 



where are plaquettes entirely in space and are plaquettes in one space and one 
time direction. The equations of motion in the temporal gauge follow immediately; 

Ui{x,t + At) = Ei{x,t + At/2)Ui{x,t) , (23) 
^Tr - iT^E,{x, t + At/2) = ^Tr - ir^Eiix, t - At/2) - {At)'^DF^{x) . (24) 

The only differences between the two update algorithms are that matrix exponenti- 
ation is replaced by inverting (l/2)Tr — zr", and the definition of energy is changed 



from E''E°'/2 to (1 - ^ {At)"^ E°- E'') / {Atf . These changes are quadratic in \E''\At, 
which is no larger than the natural level of error of the leapfrog algorithm, so the 
difference is unimportant. 

For finite At there is no identically preserved quantity which can be understood 
as a system energy; but a close reasonable approximation is the quantity 

O ^ X,i ^ 

for the first update algorithm and 
Energy (t) = ^l-iTrf/n(t) 

+ E (l - ^TrE,(a;, t + At/2) + 1 - ^TrE,(x, t - At/2))(26) 

for the second. In each algorithm, the energy is preserved in the mean, with fiuc- 
tuations suppressed by (At)^ and further washed out because they generically add 
incoherently between different degrees of freedom. The central value is absolutely 
stable, because the update algorithm has two important properties: (1) it bases the 
evolution of the system on the minimum necessary amount of state information; and 
(2) it is exactly time reversal invariant (to evolve the system backwards, flip the sign 
of E in the flrst version or take Hermitian conjugate in the second). It is therefore 
impossible that the energy should tend to rise exponentially, since under time rever- 
sal it would tend to shrink exponentially, but the algorithm is the same, further, by 
good fortune the roundoff errors which are inevitable in a numerical implementation 
add incoherently and do not tend to increase the energy. The stability of the system 
energy is necessary for long Hamiltonian evolutions and will also allow checks that 
the chemical potential method is behaving properly. 
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2.3 Magnetic field and Ncs 



d'x ' ■ (27) 



To track the evolution of Ncs in the realtime evolution of the Yang-Mills fields de- 
scribed in the last section it is necessary to develop a definition of dNcs/dt. (It is 
impossible to define Nqs directly on a lattice, see 0). In the continuum the relation 
is 

dt J ^ 87r2 

We have already found lattice analogs of Ei and ~ -Bj. Ei is associated with a 
link, and Fj^ is associated with a plaquette in the j, k directions (cyclic permutation 
implied), and the simplest definition of dNcs/dt which is invariant under the cubic 
point group is 

where 8n means a sum over the 8 plaquettes which are orthogonal to the link x,i 
and touch one of its endpoints. The plaquettes should be traversed according to the 
righthand rule, starting and ending at the point along the link, see Figure and 
parallel transport of the plaquettes which touch the endpoint of the link back to the 
basepoint is implied. In the discrete time, leapfrog algorithm evolution of the system 
the appropriate choice for ANcs from one update is 

x,i 

where E°'{x) should be replaced by (l/2At)Tr — iT°-Ei{x) for the algorithm used by 
Ambjorn et. al. 

The coefficient is 2??^ rather than 87?^, as it is in the continuum, because the 
continuum equations are based on fields defined in terms of r°/2, and throughout I 
use t"" in the discrete system. Similarly the g^, which can appear in the continuum 
equations depending on the normalization of A^, does not appear here, but will appear 
in the definition of lattice temperature instead. 

The definition of ANcs above corresponds to that used in P, 0. The definition 
used in was simpler, but not invariant under the cubic point group; so although it 
does correspond to Ncs iii the continuum limit it is contaminated by nonrenormal- 
izeable operators already at dimension 5, whereas the symmetric definition is only 
contaminated at dimension 6. 



2.4 Chemical potential for Ncs 

A particularly efficient way to study the rate of topology change in the symmetric 
electroweak phase is to apply a chemical potential for Ncs and to measure NcsT / [i. 
To do so it is necessary to develop a definition for the chemical potential, since there 
is no good lattice analog to the continuum definition of Ncs- The most naive solution 
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is that a chemical potential modifies the dynamics to 



= ,Enx)T'U,{x) (30) 

^ . -DFti.)-^-^ ,31) 

which look just like the corresponding continuum equations. However, this does 
not work, because the B term in the second equation does not preserve the Gauss 
constraint, ie 

J2iBtix)-Btix-^))^0 (32) 

i 

(where Bf{x — i) should be parallel transported as an adjoint object along the link 
leading to x). The corresponding continuum expression, DiBi, is zero, by the Bianchi 
identity; but the proof of the Bianchi identity depends on continuity of the A fields, 
and the notion of continuity is lost on the lattice. Eq. (|3^) is zero in abelian gauge 
theory. For SU(2) in the low temperature or small a limit, the violations are higher 
order in temperature than B, but nevertheless the failure of Eq. (|3^ ) to equal zero 



means that the proposed implementation will spoil the Gauss constraints, exciting 
unphysical modes. As I discuss in 0, this makes the results unreliable. 

The solution is to take very literally the constraints. What they mean is that we 
are not really living in the phase space of all E, U, but on the constraint subspace, 
where the E satisfy Gauss' law. The only actual dynamical fields are the linear 
combinations of E which are orthogonal to the Gauss constraints; if we could, we 
would express the system and its evolution in terms only of these. 

Consider the most obvious basis of electric fields, where a basis element is a 
specification of a lattice point x, a direction i, and a Lie algebra element a. 1 will 
write this basis as {Ea}- This basis has the important property that, in terms of the 
inner product defined by 2He, it is orthonormal. I write the Gauss constraints in 
this basis as G„ = CapEp = 0, where note that the first index of c only runs over the 
3A^^ Gauss constraints [N the number of lattice points on a side of a cubic lattice), 
while the second runs over the 9N^ independent electric field basis elements. (Sums 
on repeated greek indicies are always implied.) Note also that some of the coefficients 
Cai3 are U dependent, because the E fields on the links leading into a point must 
be parallel transported to the point. Now assuming that all Gauss constraints are 
independent^ we can find a different orthonormal basis in which Cq,^ is nonzero only 
for P < 3N^, so the basis can be partitioned into two subsets, the E'^, which are an 
orthonormal basis for the Gauss constraints, and the E*, which are the dynamical 
degrees of freedom. Now the definition of the discrete dNcs/dt can be written as 

2ti^'^^^ = E' ■ B + E* ■ B . (33) 
dt 



^in the abelian theory with periodic boundary conditions this is in fact not the case; there 
is precisely one linear dependence between the Gauss constraints. But for generic U the Gauss 
constraints in the nonabelian theory are independent. 
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or equivalently 



Since we are only interested in the situation in which all E"^ are exactly zero, we may 
as well define dNcs/dt as 

2n'^ = E*.B, (34) 

and hence the action of the chemical potential should be to modify the E* in propor- 
tion to the parallel component of B, but to leave the E*^ alone. In other words, if we 
write the orthogonal transformation between the E^, E* basis and the Ea b 

E^ = apaEa , Ep = bfSaEa (35) 

where orthogonality ensures 

6a-y , (36) 

then the update rule should be 

= —DFa — -^^apaap-yB^ (37) 

— = —DFa — (Ba — bpabp-yB^) . (38) 
at Zti^ 

This rule identically preserves the Gauss constraint while refiecting the meaning we 
want for a chemical potential for NcS) that it changes electric fields in proportion 
to the parallel magnetic field component. Particularly important, it will satisfy the 
property that the change in the system energy 

6H = E ■ 6E = E* ■ 6E* = - ^^ (39) 

is —fiNcs, as it should be. This allows a check of the algorithm by comparing 
the accumulated energy to Ncs- (Note that this is where it is essential that the 
algorithm without chemical potential should preserve system energy, in the mean, to 
high precision.) 

The problem with this update rule lies in finding the orthogonal transformation to 
the basis which splits into constraints and dynamical fields, which requires diagonal- 
izing the matrix Cap- Even if the diagonalization procedure could be done efficiently, 
using it via Eq. (|37| ) or (|38|) would take 0{N^) computations, which is prohibitive. 
However, notice that Eq. ( PB]) shows that the correct algorithm is to perform the 
naive implementation, but then to orthogonally project back to the constraint sur- 
face. This can be done approximately, but very accurately, by an algorithm with 
0{N^) computations per update. 

The idea is the following the most general way of moving orthogonally to the 
constraint surface is to change Ea by 

6Ea = npcpa , (40) 

with K, completely arbitrary. (It is essential that the motion be strictly orthogonal to 
the constraint surface so the dynamical fields are not changed.) The right choice of 
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K is the one which will change E so as to satisfy all Gauss constraints. A reasonable 
first guess, which is guaranteed to reduce the sum of squares of Gauss constraint 
violations, is = —•yCaisEp, < 7 < 1/6, which is a step in a dissipative or 
quenching algorithm for the Gauss constraints. (The condition 7 < 1/6 is demanded 
by stability.) A crude approximate projection algorithm is to do one application 
of this quench at each time step. A more refined algorithm is to guess that the 
correction should be almost the same as in the previous timestep, and to first apply 
= mK/3,previousstep, where m < 1 ensures algorithm stability and K^^previous step is 
the sum of all corrections applied at the previous timestep. Then the dissipative 
step is applied twice, with different stepsizes 71 and 72. I find that using the choices 
m = 1 — 1.6At, 7i = 5/24, 72 = 5/48 is stable and highly efficient. (Note that 
when two applications of the quench algorithm are performed in a row, the stability 
criterion changes, see [Q.) Applying this approximate projection algorithm requires 
about as many floating point operations as determining the DF^, so it is not a great 
strain on algorithm efficiency. It is also close enough to an exact projection that the 
residual failure has virtually no effect on the results 0. 

2.5 Thermalization 

Finally I discuss the thermalization of the lattice Yang-Mills system developed in 
this section. What a thermalization algorithm needs to do is to draw a set of E, U 
from the constraint submanifold of the phase space of E, U, choosing with weight 
exp —jSiH times the natural measure on the constraint surface. This definition of 
I3l corresponds to (3l = 4/9continuum/(5'^a), with a the lattice spacing in the same 
physical units as /^continuum- (The factor of is the usual Yang-Mills wave function 
normalization and the factor of 4 arises because I work in terms of r° rather than r'*/2, 
which are used in the definition of the continuum fields.) It is possible to perform the 
thermalization for a general constraint manifold, using a set of Langevin equations 
0; but by using the special properties that the subspace of E fields at a fixed value 
of f/ is a vector space, that the constraints are linear, and that the Hamiltonain is 
separable into a part which depends only on U and a part which is quadratic in E at 
fixed U, it is possible to find a simple and efficient thermalization algorithm 0. 

The idea is as follows; while the U dependence of the Hamiltonian, Hu, is quite 
nontrivial, all of its complexities appear in the equations of motion; and the equa- 
tions of motion preserve the thermal probability distribution, while mixing excitation 
between the E and U fields. Since we must implement the equations of motion any- 
way, it is best to use them to thermalize the system. All we need is a thermalization 
algorithm for the E fields at fixed t/; we then repeatedly draw the E fields from the 
fixed U thermal distribution, evolve the system to mix the thermalization between E 
and U fields, and throw away the E fields and repeat. The idea is identical in spirit to 
the "molecular dynamics" Monte-Carlo algorithm of Euclidean lattice gauge theory. 
The only complication here is that the E fields must satisfy constraints. 

The thermalization of the E fields is possible because the E fields at fixed U form 
a vector space, the constraints are linear, and He is quadratic in E. In terms of 
the orthogonal basis which partitions into E* and E'^, the Hamiltonian is just He = 
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+ we just need to choose each £"* from the Gaussian distribution 

Vm = V^exp -(3L{E:r/2 (41) 

and set each E'^ to zero. The easiest way to do this, given the difficulty of finding 
the basis which partitions into E* and E'^, is to set both equal to the Gaussian 
distribution, and then orthogonally project to the constraint surface. To get the initial 
Gaussian distribution, I simply draw the E^, from the same Gaussian distribution; 
because the desired basis is gotten by orthogonal transformation from this natural 
basis, the resulting distribution will also be Gaussian in the basis of E* and E'^. The 
orthogonal projection is performed by using the quenching step of the algorithm of 
the last section repeatedly until the projection is complete. I find that order 100 
applications reduce the Gauss constraint violation to the level of roundoff error (for 
single precision numerics). 

In practice I find this algorithm to be very efficient. For instance, the total energy 
of the system approaches the thermal average value by a factor of 1/2 on each E 
randomization, if the Hamiltonian evolutions are taken to be of length ~ (3l] within 
7 or 8 applications the energy has reached the level of random fiuctuation expected 
in a system of finite size, and the values of Wilson loops equal their thermal averages 
to within expected thermal fiuctuation. Two or three applications are enough to 
re-randomize a thermalized distribution. 

This supplies all necessary equipment for the numerical study of classical 3-D 
Yang-Mills theory under the minimal action. 



3 Improved Hamiltonian 

In this section I construct an improved Hamiltonian for the classical lattice Yang- 
Mills theory. I also present an improved definition of dNcs/dt and discuss how to 
implement a chemical potential for this quantity, and how to thermalize the system. 



3.1 a false start 

The simplest guess for an improved Hamiltonian takes inspiration from the example 
of scalar theory presented in the introduction and replaces the gradient term, Hu, 
with an improved term. 



Hu„e^ = o E 1 - oTrf/n - T7T E 1 - oTrf/^ , (42) 





which is known from the study of Euclidean lattice gauge theory Here J2nD refers 
to a sum over all 1 x 2 rectangular plaquettes, of which there are 6A^^. The evolution 
of the system is the same as in the previous section, except that the definition of 
DF^{x) changes to 

OfTi^) = I (E ^Tr - - 1 ( E - (43) 
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where 12 CZH refers to the 12 rectangular plaquettes containing the hnk Ui{x), tra- 
versed so as to begin and end at x and to contain Ui{x) and not Uj{x). 

To see that this proposal fails, it is useful to consider a sinusoidal excitation in 
a single Lie algebra element, and to compute the dispersion relation. Consider the 
infinite volume limit of the lattice and start from the initial conditions 

Ui{x) - / , Ei{x) = er'^Ei sin (^k ■ {x + ^i)) (44) 

where e is infinitesimal. (The coordinate in the argument of sin is the center of the 
hnk.) The energy per site of this arrangement is J2i{Ei)'^e^ /2. An infinitesimal time 
St later some of this energy has appeared as magnetic energy; the fraction is St'^cu'^, 
and cu gives the dispersion relation. At time aSt, the U fields are 

Ui{x) I + ir''Ei{x)6t (45) 

and a lengthy but straightforward calculation gives the energy from square plaquettes 

as ^ ^ 

({E,S2 - E2S,f + {Eiss - Essif + (E3S2 - ^2^3)') , (46) 

(where I introduce the shortand notation Sj = 2 sin(a/ci/2)), the energy from rectan- 
gular plaquettes as 

((8 -si- sl){E,S2 - E2S,f + (1 ^ 3) + (2 ^ 3)) , (47) 
and the Gauss constraint as 



s^E^ + S2E2 + s^E^ = . (48) 

by adding the square of the Gauss constraint, which is zero, to the energy in magnetic 
fields, I readily find that the dispersion relation is 

aW = {sl + sl + sl + ^{s\ + s\ + st)) (49) 

+ ^2 ^ ^ El (]^^?(52 + ^3 - si)) + (1 ^ 2) + (1 ^ 3) (50) 

+ E^ + m + Ei {^'^'^^'' + ^^)) + (1 - 3) + (2 ^ 3) . (51) 

The contribution to uj^ from the first line is A;? plus zero times terms of form a^k^, 
plus 0{a'^k^), and is therefore the correct dispersion relation, up to corrections which 
vanish as in comparison to the leading term. The other terms, which vanish when 
EiSi — E2S2 = E3S3 = 0, spoil the dispersion relations by introducing polarization 
dependent 0{a'^k'^) contibutions, so the evolution is not "improved". 

An even simpler way of seeing that the Hamiltonian is not "improved" is to con- 
sider the system's thermodynamics, which are determined by integration of the fields 
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over the weight exp —jSiH. Here it is convenient to follow [§ and introduce Lagrange 
multipliers ^o(^) Gauss constraints, so the Hamiltonian is 

H = Y. [^^^^^^^^ + E ^^S(^) (E Et{x) - Et{x - + /f^,^p.oved (52) 

where again parallel transport of E^{x — i) to the point x is to be performed along 



the link leading into x. In terms of the notation of Sec. 2.4 I can write this as 



,improved • 

(53) 

If we are only interested in the values of gauge invariant operators constructed from 
the link matricies then we can perform the integral over the E fields, which is now 
Gaussian, and gives 

jimproved • 

(54) 

Note that AQpCpaEa is (performing the sum on f3 and expressing the sum on a) 
J2x,iEi{^)i^oi^) "^0(^ + 0)) parallel transport of the forward Aq implied, as usual. 
So, the Aq expression above can be untangled to read 

H = l E(^0(a;) - A^oi^ + m^Ux) - + t)) + HuMpro.ed , (55) 

x,i 

which is the lattice action of the finite temperature quantum theory in the approx- 
imation of dimensional reduction, with zero bare Debye mass. It has the improved 
plaquette action, but the unimproved scalar (Aq) gradient term, as a comparison 
with the scalar theory presented in the introduction will show. Hence the Hamilto- 
nian has not been "improved," even at the level of thermodynamics, and this failure 
has nothing to do with the choice for an improved magnetic term. 

3.2 Hamiltonian and Gauss constraints 

The problem with the approach above is that, while I have removed nonrenormalize- 
able operators such as a?'Yl,ij{DiFijY{DiFijY , I have left in the nonrenormalizeable 
operator ^^{DiEiY{DiEiY . An easy way to see why is to recall how the minimal 
Hamiltonian system could be interpreted as arising from an action on a lattice in both 
space and time, consisting of all 1 x 1 plaquettes, in space and time. The E"^ term in 
H is just the sum over the 1x1 plaquettes with time links. In the previous section, 
I included 1x2 plaquettes only if each direction were in space, but continued to use 
a kinetic E"^ term derived only from the square plaquettes. Instead I should include 
rectangular plaquettes here too. The result is that the continuum time Hamiltonian 
should be 

_ 5 ^ {Et{x)f 1 ^ {2Et{x)f 1 ^ (Ef (x) + Efjx + i)f 

J^E.irapYoved — g ^ ^ 2 12 2 ' 

i^x i,x i,x 
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where the first term comes from square plaquettes, the second comes from plaquettes 
which stretch 2 units in time, and the third comes from plaquettes which stretch 2 
units in space, and I have aheady taken the small time step limit. As usual, E'^{x + i) 
must be parallel transported along the straight path. The second term should be 
absorbed in the first, changing its coefficient to 4/3. I have retained my earlier 
definition of E°-{x), so in the Aq = gauge the update rule for the U is 

^^ = ir^Et{xmx). (57) 
The Hamiltonian is conveniently summarized by defining the operator M, 

{METdx) = ^Etix) - ^{e^{x - z) + E^{x + t)) (58) 
(parallel transports again implied), in terms of which the Hamiltonian is 

H = Hu ,improved 

+ -E^M^pEp . (59) 

The Gauss constraint is derived as before, by varying Aq{x)\ it is the sum of all 
plaquettes which traverse the time direction and terminate at x, with appropriate 
signs: 

^ ' {Et{x) + Et{x + i)) - (Etix - z) + Efix - 2i)) ) ) , (60) 
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or in somewhat more compact notation, 

CpaMayE^ = . (61) 

Now I will address the question of whether this Hamiltonian is improved in the 
sense of replicating the continuum evolution for smooth fields without contamination 
from dimension 6 operators. There are two general types of dimension 6, nonrenor- 
malizeable operators which can appear in cubic point group invariant lattice imple- 
mentations of Yang-Mills theory; there are those with two powers of the field strength 
tensor, such as J2ij{DiFij)"'{DiFij)°- , and terms which contain three powers of the field 
strength, like fabcF^jF^j^F^^. In the latter term, the indicies of the F must differ, so 
that they represent curvatures in different planes, since fabcF^jFj'j automatically van- 
ishes. Since all the plaquettes used in Eq. (^) are planar, terms involving three 
powers of the field strength will not appear at dimension 6. I can therefore restrict 
attention to the other sort of term, which appear already in the Abelian theory and 
which modify the dispersion relations of weak disturbances which vary sinusoidally 
in space. It is sufficient, then, to check whether the dispersion relation for uj"^ is free 
of /c'^a^ terms in the small k expansion. It is again sufficient to start a sinusoidal 
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excitation with U = I everywhere and see how quickly the energy is transferred into 
the magnetic term of the Hamiltonian. For 



Ei{x) = eT''^;^ sin(A; • {x + i/2)) (62) 
the electric energy is times 

\ {El + El + El) - ^ ((4 - sl)El + (4 - sl)El + (4 - sl)El) , (63) 

with the first term from the single E term in the Hamiltonian and the second term 
from the two E term in the Hamiltonian. Again I use the notation Sj = 2sin/ci/2. 
The new Gauss constraint is 

(l + ^5?) E^s^ + + ^si) E2S2 + (1 + ^stj Esss = , (64) 

and the magnetic energy at time aSt is still 

(1 + ^(^1 + ^2)) (^1^2 - E2S,f + (1 ^ 3) + (2 ^ 3) (65) 

as in the previous subsection. Squaring the Gauss constraint, which equals zero, and 
adding it to the magnetic energy term, I find that the oscillation frequency is 

aW = Us^ + ^st) (66) 

^ 1 slsl{E,S2 - E^sif + (1 ^ 3) + (2 ^ 3) 
144(i + £|)^2+(i + ||)^|+(l + £|)^|- 

The first term here is 

8 1 

- ^(1 — cos a/cj) ^(1 — cos2aA;i) (68) 

i i 

which is the usual dispersion relation for an improved, scalar Laplacian. It contributes 
no a?k^ terms to o;^, as can be readily verified by expanding the cosines as power series. 
The second term is polarization dependent, but has 6 powers of 2smak/2, so it only 
starts to contribute at order a'^k^, and the dispersion relation is therefore free of the 
unwanted a'^k^ terms, and the Hamiltonian is free of dimension 6 nonrenormalizeable 
operators. Also note that the new Gauss constraint demands that E be transverse to 
a higher order in k than the unimproved Gauss constraint does. 

I can also check the thermodynamics as before. Introducing Aq{x) as a Lagrange 
multiplier for the constraint, the thermodynamics is governed by the probability dis- 
tribution exp —(3lH , with 

H = ifi7,improved + Ij^EaMafjEp + iAQpCp^M^aEa (69) 
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which, on performing the Gaussian integration over E, becomes 

-^(7,improved + ^ApCpaMa-yCs^Aos . (70) 

The Ao term untangles to 

^{Ao{x) - Ao{x + i)) [^Ao{x + 2i) - ^Ao{x + i) + ^Ao{x) - ^Ao{x - i)) , (71) 

parallel transports again implicit. The derivative term for Aq is precisely the improved 
derivative term for the scalar theory discussed in the introduction, made covariant. 
Hence the Hamiltonian does give improved thermodynamics. This also points out a 
relation between the improvement scheme for the plaquette action and for the discrete 
Laplacian. 



3.3 Numerical evolution of the improved system 

Rather than derive continuum time equations of motion from the Hamiltonian, I will 
skip to finding a discrete time evolution algorithm. Naively, I should base this on the 
action 



S 

Se 



Se — Su 



Su - E 



x.t 



^El - ^Tr^-. - ^El - - ^El - ^Tr^-0. 



^ Vl - ^Trt/p.. - - Vl - -TrC/^.. - - Vl - ^Trt/a. 
3 ^ 2 12 -f^ 2 12 2 ' 

i>] i>] i>j 



(72) 



which includes plaquettes which go two steps in time. These plaquettes prevent 
a nonrenormalizeable operator of form (At)^(DoFoi)"(-Do-^oi)" and should make the 
evolution fourth order in stepsize. This is unnecessary if we take At « 1, and is also 
undesirable since the inclusion of plaquettes which stretch two steps in time makes the 
numerical evolution unstable. The reason is that the update rule derived by varying 
one link then determines U based on the 4 previous time slices, rather than just 2 
as the minimal algorithm did; since the Yang-Mills equations are second order, this 
is too much state information, which opens the possibility (which is realized) of an 
unstable growing mode. Instead I choose At small enough to make those plaquettes 
unnecessary, and choose 



E 

x,t 



3 2 °^ 12 ^ 2 ° 



(73) 



As before, I will define Ei{x,t + At/2) as the square plaquette, so in the gauge 
Uo{x) = /, 

Ui{x, t + At) = Ei{x, t + At/2)U^{x, t) . (74) 
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The Gauss law is that the sum of Lie algebra traces of all plaquettes containing a 
time link be zero, which in the same gauge is 



= E[^QTr-2r'^f/,(x,t + At)f/;(a;,t)- 



(75) 



1 /I 



12 V2 



^Tr - ir'^Ujix - i, t)Ui{x -i,t + At)^ - 
-Tr - ir^Uiix, t + /\t)U^{x + i,t + At)Uj{x + i, t)Uj{x, t)- 



1. 



Tr - ir^'Ujix - i, t)Uj{x - 2i, t)Ui{x -2i,t + At)U^{x -i,t + At) 



or 



E 



4 /I 



3 V2 



-Tt - iT''Ei{x,t + At/2)- 



-Tt - ir^Ulix - i, t)Ei{x -i,t + At/2)Ui{x - i, t)) - 



1 n 



12 V2 



-Tr - ir^Eiix, t + At/2)Ui{x, t)E,{x + i,t + At/2)uj{x, t) - 



^Tr - tT^Uj{x - i, t)Ul{x - 2i, t)Ei{x -2i,t + At/2) x 
Ui{x - 2i, t)Ei{x -i,t + At/2)Ui{x - i, t)^ 



(76) 



Implicitly parallel transporting using the connection at time t, this is 

= E 



4/1 1 

- f -Tr - ir^'Eiix, t + At/2) - -Tr - ir"^,(a; - t + At/2) ) - 



1 /I 



12 V2 



-Tr - ir^Eiix, t + At/2)Ei{x + i,t + At/2)- 



^Tr - iT^Eii^x -2i,t + At/2)Ei{x -i,t + At/2)^ 



(77) 



whereas if the parallel transports are performed using the connection at time t + At, 
this is 

= E ^(^^Tr-iT''Ei{x,t + At/2)-^Ti-iT''E,{x-i,t + At/2)^~ 



1 /I 



12 V2 



-Tr - ir"Ei(x + i,t + At/2)Ei{x, t + At/2)- 
-Tr - iT''Ei{x -i,t + At/2)Ei{x -2i,t + At/2) 



(78) 



The equations of motion, derived by varying one spatial link, are also somewhat 
more complicated. In terms of the definition of DF^[x) appearing in Eq. (P3|), 
and implicitly parallel transporting objects using the connections at time t and the 
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straight line paths, the update rule is 

lQTr-zr'^E,(x,t + At/2)) - 
^ Qxr - iT^Ei{x, t + At/2)Ei{x + i,t + At/2)- 

^Tr - iT''Ei{x -i,t + At/2)Ei{x, t + At/2) 
1 {h:i-iT''E,{x,t- At/2) 
^ Qxr - iT^EiiyX + i,t- At/2)Ei{x, t - At/2) + 

^Tr - ir^'Eiix, t - At/2)Ei{x - i,t - At/2)^ - {AtfDF^{x, t) . (79) 

Writing this equation in the convenient shorthand 

{ME)1{x, t + At/2) = {ME)1{x, t - At/2) - {AtfDF^{x) , (80) 
one can see that the Gauss constraint is 

J2i{ME)^{x, t + At/2) - {ME)^{x -i,t + At/2)) = (81) 

i 

and that the difference between the value of this quantity at times t + At/2 and 
t- At/2 will be 

- (At)' j:{DFt{x, t) - DFtix - t, t)) (82) 

i 

which vanishes identically, because each plaquette which appears in the evaluation 
of one term appears with opposite sign in the evaluation of another term, as before. 
Hence the Gauss constraint will be preserved identically by the evolution. 

Eq. (^) does not give the Lie algebra content of the electric fields, (l/2)Tr — 
iT°'Ei{x), directly, but gives the result of a nondiagonal, nonlinear transformation 
on them; to implement the update it is necessary to invert this transformation. The 
inversion problem is not too bad (for At -C 1), because the expression on the lefthand 
side of Eq. ([79|) is close to diagonal in the E fields and the inversion can be performed 
iteratively, ie 

^ Qxr - iT^'Eiix)^ = (r. h. s. of Eq.(79)) + (83) 

^ Qxr - iT''E,{x)Ei{x + t) + ^Tr - ir'^E.ix - i)E,{x) - 2^Tr - ir'^Eiix) 

can be solved by making a guess for E, substituting into the righthand side, and 
getting a better guess for E. I have written the formula so that the righthand side 
depends almost only on E at the points 1 space ahead and behind of the point of 
interest, so the convergence rate is doubled by alternately updating even and odd 
lattice points. 
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I have implemented this algorithm numerically. At single precision, choosing the 
initial guess for Ei{x,t + At/2) to be Ei{x,t — At/2), it is only necessary to iterate 
over the lattice 4 times before the remaining correction comes on order roundoff error. 
If the inversion were conducted completely to infinite precision then the algorithm 
would not gain or lose energy in the mean, because the update rule is time reversal 
invariant and only depends on the minimal amount of state information. For a single 
precision implementation, because the iteration converges from above, the roundoff 
errors systematically boost the energy, so that it rises by about 2 parts in 10^ per 
time step; but this level is small enough to be of no concern if the total length of 
a simulation is kept under t ~ 5000. By going to double precision and making the 
timestep very small to reduce the fluctuations in system energy I was able to verify 
that the algorithm does preserve energy in the mean. I also observe that violation of 
the Gauss constraint is only generated by roundoff error, as in the evolution with the 
unimproved action. As a final check I verified that the implementation produces the 
expected dispersion relations when the initial state is a weak sinusoidal excitation in 
one Lie algebra direction. 

3.4 Magnetic field and A^^^ 

The next step in the program is to find the appropriate definition of dNcs/dt for 
the improved Hamiltonian. The approach is similar to the unimproved case, where 
J{dNcs/dt)dt is defined as 

where the first sum is over all points in the 4 dimensional grid, the second is an 
average over the 4 plaquettes in the jk plane {ijk a cyclic permutation) which touch 
the point, and the last is a sum over the 4 plaquettes which proceed one step forward 
or backward in the time direction and a step forward or backward in the i direction 
from the point. That is, it is an average over the electric fields leading into and out of 
the point, taken a half timestep before and after the spacetime point. In the previous 
section I rearranged this sum into a sum over the electric fields of a contribution 
coming from each end of the link on which the electric field resides, which is how the 
sum over 8 square plaquettes arose. 

To modify this definition so as to remove dimension 6 contamination, I need to 
replace the sum over square plaquettes with a sum on square and 1x2 rectangular 
plaquettes. Again, because the plaquettes are planar and the improvement need only 
correct dimension 6 errors (dimension 4 in the electric or magnetic field), terms cubic 
in the field strength will not arise and 1 need only consider weak sinusoidal fields in 
one Lie algebra direction. In a small k expansion, the choice of spatial plaquettes 
should not contain any k^ term. 

The contribution from the 4 square plaquettes around a point x and in the 1, 2 
plane for Ui{x) = 1 + eir" sm{k ■ x + ki/2) is 

AciC2{E2Si - EiS2)e'^ sm{k ■ x) , (85) 
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where q = cos/cj/2. The contribution from the 4 rectangular plaquettes which have 
midpoints at x is twice this, and the choice of whether to use square plaquettes or mid- 
points of rectangular ones is ours (though it may become forced when improvement 
is carried up to higher order). The contribution from the 8 rectangular plaquettes 
which have a corner at x is 

4ciC2(4 -si- si) {E2S1 - EiS2)e^ sin{k ■ x) . (86) 

The correct combination should produce no k'^ term. The requirement is solved by 
using 5/3 the sum over square plaquettes (or 5/6 the sum over rectangular ones with 
midpoints at x) plus —1/6 the sum over rectangular plaquettes with corners at x. For 
the electric field it is not necessary to use rectangular plaquettes which stretch two 
units in the time direction, because the field variation in this direction is suppressed 
by a factor of (Af)^, and the right combination at the point x,t (making parallel 
transports using the links at time t) is 

:T''Ei{x) + ^Tr - iT''Ei{x - i)\ {t + At/2) - (87) 
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This differs only by terms of order (At)'^ from the same expression with the trace 
taken over each E in the two E terms, rather than over the product. This is because 
the order At correction arises from a commutator term which enters with opposite 
sign in the t — At/2 and t + At/2 contributions. Since there are already inevitable 
0((At)^) errors in the evolution of the system, I am free to choose the more convenient 
definition without worrying about the order (At)^ change it will make in the results. 
Choosing to trace each E separately, I can rearrange the definition to be 

1 /I. 



2^ANcs = J2Bt{x,t)-(^-Tr-iT''Ei{x,t + At/2)+ 



lTr-zr"E,(a:,t- At/2)) 
Bf{x) = zlF^^{x + 2i) + ^F;^,{x + i) + ^F;^,{x)-^F^,{x-i) (89) 
^m^) ^ I E - i^^Uu,, - ^ E ^Tr - zr^U^^, (90) 

where all definitions should be clear from earlier usage and parallel transports of Fj^. 
are implied along the straight line paths to the point x. 

This gives a usable definition for the evolution of Ncs which is free of dimension 
6 contamination. 
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3.5 Chemical potential for Ncs 

Now I should apply a chemical potential for Nqs in order to study its motion due 
to infrared physics. In analogy with the unimproved Hamiltonian case, the chemical 
potential term should modify the evolution of the E fields in a way which makes the 
overall energy change by —fiNcs- In the unimproved case this was accomplished by 
making 

so the change in the electric energy was 



SE^ = -^B^, (91) 



EJE^ = -^^^1^ = -^6Ncs . (92) 

Then I modified the definition so that the linear combination of E fields representing 
a Gauss constraint violation were not excited after all, by orthogonally projecting 
to the constraint surface. This did not change the energy because these modes had 
started out zero anyway. 

In the improved case, it is easiest to begin by thinking about the continuum time 
system. The change in Ncs is 

EaBadt 

dNcs = „ 2 ' ^3 
27r^ 

but the shift in the energy is SEaMapEis, so 6E should be 

before worrying about preserving the Gauss constraint. Note that the evolution 
equation for E already includes a term like —M~i^ DF/^dt, so the magnetic field appears 
in the same place as DF^{x). This describes a possible implementation in the discrete 
case as well, namely that Eq. (|79|) is modified by replacing —DF°'{x) by —DF°'{x) — 
/i5f(x)/27r2. 

Another way of arriving at this result is to note that the Hamiltonian defines a 
metric on the space of E, and that in the improved Hamiltonian the natural basis E^ 
is no longer orthogonal. Rather, the orthogonal basis is (M^/^)q,^£^^. In terms of this 
basis, the definition of Ncs is 



dNcs _ m-')apEp){M-2)^^B^ 



dt 27r2 
and so the change in the electric fields should be 



(95) 



S{{M2)^f^Ep) = , (96) 



which leads to Eq. (B 
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Now although the new definition of the magnetic field should be free of the leading 
order of operators responsible for it not satisfying the Bianchi identity, so the satisfac- 
tion of the Gauss constraint should be better in the infrared, the ultraviolet violation 
of the Gauss law is just as bad as in the unimproved case, and it is again necessary 
to change the definition of the action of a chemical potential to prevent the violation 
of the Gauss constraint. In analogy with the unimproved case, the correct modifi- 
cation is to orthogonally project back to the constraint surface. Again, orthogonal 
should mean with respect to the metric defined by He-, so the most general motion 
orthogonal to the constraint surface is 

5Eo, = l^pCfSa not SEa = KpC0.yM^a ■ (97) 

For small k, the former changes the system energy by 

6He ^ E^M^p5Ep = E^M^pc^pKp , (98) 

which is suppressed by a power of the adherence to the Gauss constraint, whereas the 
latter changes the system energy by 

E^Mo^pMp.c^.Kp , (99) 

which is nonzero even if the Gauss constraint is satisfied; so the modification to E is 
not orthogonal to the constraint surface. 

A nice way to understand this conclusion is to remember that M^^'^E are orthog- 
onal; the Gauss constraint is 

= {c^p{m"^)pMm'^UEs) = 0, (100) 

so the most general correction to the electric field which is orthogonal to the constraint 
surface should be 

5{{M^^pEp) = K^{M^)aeC^, or 5Ea = Kf^CfSa (101) 

which is the same as Eq. (p7|) . 

I can now find an approximate projection algorithm in strict analogy with the 
unimproved Hamiltonian case. A step in a dissipative or quench algorithm for the 
orthogonal projection is 

6Ea = KpCfia , tifS = -'^cp^M^^E^ (102) 

where stability now demands 7 < 1/8. This step is orthogonal to the constraint 
surface, and repeated application is guaranteed to be a projection; a single application 
at each step is enough to keep the departure from the constraint surface under control. 
A better algorithm, again in analogy with the unimproved system, is to make an initial 
correction with k equal to m < 1 times the total k of the previous step, and then to 
apply the dissipative algorithm twice, with different stepsizes. 

Note that Eq. (^) is of exactly the same form as in the unimproved case. To 
understand why, recall the origin of the Gauss constraint; it is the Lagrange equation 
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for the Aq field (or the time direction link matrix), and when it is not satisfied it means 
we have inadvertently not set Aq to zero; since the electric field only corresponds to 
the time derivative of the link matrix when Aq is zero, we should restore Aq to zero, 
and Eq. ( |101[ ) is precisely how the system changes as we modify Aq by k. 

This latter observation explains how to project to the constraint surface in the 
discrete time case. The modification of the electric fields should take the form of a 
time dependent gauge transformation at each point, with magnitude chosen to restore 
the Gauss constraint everywhere. In other words the correct orthogonal projection 
to the constraint surface is a choice of a group element g{x) at each point x which 
modifies the E matricies through 

Ei{x,t + At/2) g{x,t)Ei{x,t + At/2)Ui{x,t)g\x + i,t)U}{x,t) . (103) 

The g{x) should be chosen such that 



G^{x,t) = y: 



4/1 1 

- f -Tr - ir^Eiix, t + At/2) - -Tr - ir"E,(x -i,t + At/2) 



^ ^Tr - ir^'Eiix, t + At/2)E,{x + i,t + At/2)- 
^Tr - ir^Ei^x -2i,t + At/2)Ei{x - i,t + At/2) 



(104) 



(parallel transports at time t implied) is zero. Again the job of finding the right g is 
intractable, but a reasonable guess (the basis of the dissipative algorithm) is 

^Tr-^r"^(x,t) = -7G"(x,t), (105) 

and a better guess is to use m times the Lie algebra element used in the previous time 
step, and then to apply the above dissipative step twice with different coefficients 7. 
The stability condition is now —1/3 < (1 — c<J7i)(l —00^2) < 1 for all < a; < 16, and 
I find that a good choice is m ~ 1 — l.GAt, 71 = 5/64, 72 = 5/32. Using these values, 
the Gauss constraint is as nearly satisfied as in the unimproved Hamiltonian, and the 
relation AH = —fiANcs is satisfied to < 1% (after energy gain due to roundoff errors 
are taken into account; I checked this in a double precision implementation, where 
there is no energy gain from roundoff error) , which should assure that the systematic 
errors in the determination of the linear response to a chemical potential, due to not 
using an exact orthogonal projection, are of the same size Q. 

This provides a satisfactory implementation of chemical potential in the improved 
Hamiltonian system. 



3.6 Thermalization 

Next I extend the thermalization algorithm developed for the unimproved Hamilto- 
nian to the improved one. I will begin by showing how to thermalize the system in 
the continuum time case. The technique does not strictly work for finite At, but the 
failures are of order (At)^, and since the time evolution and the definition of Ncs 
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already introduce systematic errors of at least this size I will ignore this problem and 
find a discrete implementation which is correct up to these stepsize ambiguities. 
In the continuum time case, I want to thermalize the system with Hamiltonian 



given in Eq. (|59D, subject to the Gauss constraints given in Eq. (^). The subspace 
of momenta at fixed U is an inner product space with 2He the inner product, and 
the Gauss constraints are linear, so I can implement the thermalization discussed in 
Subsec. p.5| , modified only to account for the new inner product; in other words I 
need only find a thermalization algorithm for the fixed U electric field part of the 
Hamiltonian and then use the equations of motion to stir the thermalization between 
electric and magnetic modes. 

If we had an orthonormal basis which partitioned into Gauss constraints and dy- 
namical degrees of freedom, we would choose the dynamical degrees of freedom to be 
Gaussian distributed and the constraints to be zero, and the electric fields would be 
thermalized. Equivalently I could Gaussian distribute both of them and then orthog- 
onally project to the constraint surface. The continuum time orthogonal projection 
algorithm is given already in the previous section, and a Gaussian distribution in 
any orthonormal basis is Gaussian in any other, so the only remaining problem is to 
find an orthonormal basis. In the previous section we saw that the basis M^^'^E is 
orthonormal. Therefore, a thermalization of the electric fields consists of choosing 
the M^/'^E Gaussian, inverting M^/^ to find E^^ and applying the orthogonal projec- 
tion by repeated application of the continuum time dissipative algorithm of the last 
section. 

The inversion can be carried out as follows: note that M is (7/6)/ plus a small 
off diagonal term, M^p = 7/6{6a/3 + fna/3); write 

K = {M--^)^p{M-2E)p (106) 
{{l{I + m))-'^U{M-^E), 

3 o 5 




map + ^Kf3 - Y^rnlp + • • •) {m'^E)p , 



in other words, expand M~^/^ in a Taylor series. The absolute magnitude of the 
largest possible eigenvalue of m is 1/7, so the Taylor series converges rather well; and 
the inefficiency of the algorithm does not matter because it is only called occasionally, 
with several (order 100) lattice updates inbetween. 

It is not completely trivial to carry this algorithm over directly to the discrete time 
setting. For one thing it is not clear that randomizing Ei{x, t+At/2) is the appropriate 
thing to do. If we did, the thermalization algorithm would not be quite time reversal 
invariant, as Tr — iT°'Ei{x, t + At/2) would then be uncorrelated with DF°-{x, t), but 
backwards evolving, we would find that Tr — iT°-Ei{x,t — At/2) was correlated with 
DF°-{x,t). Also, the E are matricies in SU(2), not Lie algebra elements, and while 
there is a mapping from SU(2) to the Lie algebra, given by l/2Tr — ir", which is one to 
one from half of SU(2) and not difficult to invert in this patch, the determinant of the 
map is not 1, ie it does not preserve the measure of SU(2). Also, the energy depends 
not on 1/2 the square of the Lie algebra element, but on 1 — l/2Tr of the group 
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element, which, in terms of Ef = l/2Tr - ir^E, is 1 - y/T^^WW ^ E''E''/2. In 
terms of the Lie algebra element, nonquadratic terms will appear in the Hamiltonian, 
both from the measure and from the way the group element is traced to get the energy. 
Therefore the thermalization algorithm I have used is, strictly speaking, wrong. 

However, I can define a Lie algebra element "electric field" at integer time steps, 
for which all of these ambiguities would appear as order (At)^ or (At)^//?^ correc- 
tions, and can in practice be ignored. (It is not even clear that there is a meaning 
for temperature beyond leading order in (At)^ in the discrete time system.) The 
definition is 

^E^ix^t) - L(Et{x-t,t) + Et{x + i,t))= (107) 



1 

At 



'^^^^^-Di^'^(x) + - f ^Tr - iT^Eiix) ) (t - At/2)- 



2 ' ' ' 3 V2 



^ ^Tr - iT''Ei{x)Ei{x - z) + ^Tr - ir^Eiix + i)Ei{x)j {t - At/2) 

(all parallel transports to be conducted with respect to the links Ui at time t), from 
which it is simple to continue the time evolution of the system. The Gauss constraints 
on either side of time t will be satisfied provided that Ef{x,t) obeys 



= J:\l{E^('^^t)-E^{x-^,t))-^(Enx + ^,t)-Et{x-2^,t) 



Ca/sMi^jE^ 



(108) 



I will thermalize this E!^ by drawing it from the thermal distribution of the Hamilto- 
nian EciMa(3Ep/2. The algorithm is precisely the continuum time algorithm presented 
above. 

To check that this thermalization algorithm is reasonable, in the sense that it 
gives thermodynamic properties which change only very weakly with stepsize when 
the stepsize is small, I have evaluated the magnetic energy and Wilson loops at 
different stepsizes. Some results for a 16'^ lattice at lattice temperature /S^, = 5 are 
presented in Table ||. For each column I have used the algorithm to thermalize an 
initial configuration and then continued to implement the algorithm for 80 E field 
randomizations, recording each observable once each E field randomization. The error 
bars are from statistics, corrected for correlations because one E field randomization 
is insufficient to fully rethermalize the system. There are clear stepsize effects which 
are consistent with the expected (At)^ behavior. The changes are significant for the 
two larger stepsizes in the table, so At cannot be taken to be this large; but the 
zero At extrapolated magnetic energy differs from the At = 0.05 value by only 0.5%, 
which will lead to a systematic error in Lyapunov exponents of order 0.5% and in 
the motion of Nqs under chemical potential of order 1.5%, which will be less than 
statistical measuring error. At = 0.05 is therefore sufficiently small that the remaining 
stepsize systematics can be neglected. Note also that these systematics will apply to 
all measurements, so they will not affect any conclusions about the existence of a fine 
lattice spacing limit. 
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quantity 


at At = 0.05 


at At = 0.20 


at At = 0.30 


Hit- ABTli?,N^) 
-'-'(7 improvedML/ j 


1 031 ± 0016 

X . W (J X _l — . W W X \J 


1 1081 ± 0015 

X . X W W X 1 . W W X 


1 2278 ± 0016 

X . ^ 1 w 1 .wwxw 


1x1 Wilsnn Innn 

X /N X VV llowli IW^LJ 


84865 + 00024 


83774 + 00021 

• Kj'^J I 1^ —1— . W W WZj X 


82087+ 00025 

.Kj I —1— .WWWiiW 


1x2 Wilsnn Innn 


72376 ± 00035 


70763 ± 00034 

. I w I wo — 1 — .wwwo^ 


68274+ 00038 


1x3 Wilsnn Innn 


6201 ± 0005 


6005 ± 0005 

.WWWtJ — 1 — .WWWtJ 


5703 + 0006 

.w 1 wo — 1 — .wwww 


1x4 Wilsnn Innn 


5316 ± 0007 

• tjoxxj _i — .www 1 


5104 ± 0007 

. 5J X Wt: —I — • W W W 1 


4768 + 0007 

.t: I W \J _l — .WWW 1 


1 X \A/il^nn loon 


4560 ± 0008 


4338 ± 0008 


3988 + 0008 


2x2 Wilsnn Innn 


5412 ± 0010 

.wt:x^ 1 .wwxw 


5222 ± 0010 

tXj^^^ 1 .wwxw 


4909 + 0010 

ij — 1 — .wwxw 


2x3 Wilsnn Innn 


4121 ± 0009 

• ^X^X 1 .WWWty 


3925 ± 0008 


3592 + 0009 

tXjxJiJ^ _l — .WWWty 


9 y 4 Wilson loon 


3149 + 0010 

.OX^<_/ —1— .WWXW 


2967+ 0009 

.^jU\Ji —I— .WWW<_/ 


2644+ 0010 

. ^ Wt:t: —I .WWXW 


9 X R \A/ilc;nn loon 


2410 + 0010 

. X W — 1 .WWXW 


2246 + 0010 

. ^ ^^W — 1 .WWXW 


1Q4Q + nnin 

• A^U^U —1— .WWXW 


.'^ v .'^ \A/i 1 c;n n In on 

'J /\ -J VV llOwli IWV-'IJ 


2843 + 001 5 

. ^ —I . W W X 


2660 + 0015 

.ijwww —1— .wwxw 


2344+ 0015 

. ^ O^t: —I .WWXW 


8x4 Wilson loon 


1975 ± 0011 

. X i _l — . W W X X 


1824 ± 0010 

. X 1^ ^j^^ 1 .wwxw 


1548 + 0010 

. X Wt: w — 1 — .WWXW 


3 x 5 AVilsnn Innn 


1377 ± 0011 

. X ' J 1 1 — 1 — . KJ KJ X X 


.1255 ± .0011 


1018+ 0010 

. X \J X \^ — 1 — . \J \J X \J 


4x4 Wilson loop 


.1251 ± .0015 


.1143 ± .0015 


.0924+ .0015 


4x5 Wilson loop 


.0848 ± .0011 


.0722 ± .0011 


.0546 + .0010 


5x5 Wilson loop 


.0467 ± .0015 


.0424 ± .0015 


.0286 + .0015 



Table 1: Stepsize dependence of thermodynamic properties in the thermalization 
algorithm for the improved Hamiltonian lattice, for a 16^ lattice at lattice inverse 
temperature I3l = 5. 

Based on this analysis, and because the evolution of the system will have stepsize 
errors of similar size, I have used At = 0.05 everywhere in what follows, unless noted 
otherwise. 



4 Numerical results 

I have implemented the improved Hamiltonian system developed in the previous sec- 
tion numerically, and used it to study the motion of Ncs under a chemical potential. 
It has become conventional in the field to quote the rate of Ncs violation in the sym- 
metric phase in terms of a dimensionless quantity k. The linear response coefficient 
to a chemical potential is 



which in lattice units is 



T{Ncs) I rp\/^ 

(/5i7r)4AiVcs 



By a detailed balance argument |[T5| , ^ this is half the rate of Ncs diffusion, 

{{Ncs{t)-Ncsm?) 



Trf = lim 



vt 



Kd{awT) , Kd = 2k^ 



(109) 



(110) 



(111) 
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number of runs 


Kjii 


3 


16 


3 


181 ± 012 


4 


16 


6 


310 +017 


5 


16 


16 


389 ± 021 


6 


16 


27 


.413 ± .024 


8 


24 


16 


.487 ± .033 


10 


30 


18 


.537 ± .036 



Table 2: Rate of motion of Nqs under a chemical potential at several lattice spacings. 



measured in 0]. 

I have evaluated the rate of Ncs motion under a chemical potential for several 
lattice spacings (lattice reciprocal temperatures), using in every case a chemical po- 
tential /i = Ij i^Li which previous work shows is easily small enough to fall in the linear 
response regime [0. I used lattice sizes such that > 3/32, in most cases, because 
it has been shown that finite lattice size effects vanish already at = 2/?^ ||^. The 
evolutions were each of length t = 400 in lattice units; for each lattice coarseness I 
used several such evolutions from independent initial configurations drawn from the 
thermal ensamble by the thermalization algorithm. The error bars are statistical; 
when there were several runs they were determined by the fiuctuations between in- 
dependent runs. 1 also computed them by assuming the statistical error comes from 
Brownian motion in Nqs, with a rate given by Eq. (|111|) ; the two estimates agree 



within error, so I used the latter estimate when there were few independent runs. 

My results are presented in Table ^ and plotted, along with some results of p, Q 
for the unimproved Hamiltonian, in Figure As the lattice becomes finer, appears 
to converge to a lattice spacing independent result. The rate of convergence is slightly 
worse for the improved Hamiltonian than for the unimproved one, so the physics which 
establishes the rate of Ncs diffusion depends on the Debye screening length being 
well shorter than the length scale of nonperturbative physics. In fact, from the data 
presented it is not entirely clear that the improved lattice results have converged to the 
fine lattice spacing limit; but the numerical demands rise as the 4'th power of jSi, and 
the improved Hamiltonian requires 5 times as many fiops to evolve the same lattice 
4-volume as the KS Hamiltonian , so I was unable to get useful data on finer lattices 
than those presented. I will assume here that the answer has converged at = 10, 
since even if the convergence rate is Debye mass limited, then the convergence should 
be as good here as in the unimproved system at ^ 8. 

For the current level of statistics the results of the improved system are compatible 
with those found using the Kogut-Susskind Hamiltonian. However, the still substan- 
tial error bars leave open the possibility that the rates differ by order 5% to 10%. To 
compare the rates to 1% would require approximately 50 times the integration time 
used here, or about 2 x 10^^ floating point operations, which is prohibitive. 

To examine whether the two Hamiltonians really do produce different infrared dy- 
namics, or whether I just have insufficiently reflned lattices and insufficient statistics, 
I have measured another quantity, the maximal Lyapunov exponent, for the improved 
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N 


improved Xmax 




improved diXmax 




2 


12 




659 ± 002 




1 318 ± 004 

X . (J X _l — •\j\J j: 




12 


4Q32 + 005 


586 + 002 


1 233 + 013 


1 465 + 005 


3 


12 


399 ± 0034 


497 ± 004 


1 197± 010 

X . X i _l — . W X w 


1 491 ± 012 

X • tJ X 1 • W X ^ 


4 


12 


2863 ± 0021 


350 ± 003 


1 145 ± 008 


1 400 ± 012 


5 


20 


2280 ± 0018 


260 ± 002 


1 140 ± 009 


1 300 ± 010 


6 


20 


1908 ± 0010 




1 145 ± 006 

X.XTItJ _l — •\J\J\J 




7 


20 


1630 ± 0027 




1.141 ± .019 




7.5 


20 




.163 ± .001 




1.222 ± .008 


8.75 


20 




.137 ± .001 




1.199 ± .009 


10 


20 




.120 ± .001 




1.20 ± .01 



Table 3: Lyapunov exponent for improved and Kogut-Susskind Hamiltonian at 
several lattice spacings. The improved data is new, the KS data is taken from [0. 
Each Hamiltonian system shows a fine lattice limit, but the value of the limit differs 
by about 6% between the Hamiltonian systems. 

Hamiltonian system. The maximal Lyapunov exponent is defined as the rate at which 
two almost identical initial configurations diverge, as measured by some gauge invari- 
ant phase space metric. It is a better quantity to use to compare the improved and 
unimproved Hamiltonian, because it depends on less equipment than the measure- 
ment of K^, requiring only the evolution and thermalization algorithms, and because 

the statistics improve more quickly, rising roughly as yJJptJpf^, whereas for the 

statistics improved as fi^Hf^NH/ {(3lt^Y/2. It will therefore be possible to compare 
the maximal Lyapunov exponents of the two theories to order 1% without excessive 
numerical demands. 

I will base my approach on that of Krasnitz 0. I draw an initial condition from 
the thermal ensamble and then perturb it slightly; then I allow the original configu- 
ration and the perturbed version to evolve under the classical equations of motion, 
periodically measuring the metric distance. At first the separation is dominated by 
transients; then there is an epoch of exponentially growing separation, and finally the 
separation saturates. It is the exponential epoch which I use. I repeat for several 
thermal initial configurations to determine and improve the statistics. 

In my implementation I start by getting a thermal initial condition by applying 
my thermalization algorithm at single precision and a stepsize of At = 0.05; I then 
convert to double precision, projecting the fields to the SU{2) manifold and quenching 
the Gauss constraint violations to the level of double precision roundoff error. I then 
copy this initial configuration and evolve it one time step, both with and without a 
tiny (10^^'^) chemical potential, to generate initial differences between configurations 
at the level of double precision roundoff error. I let the configurations evolve, tracking 
the metric distances 

De = Y: |^TrE,(x) - l-TrElix)\ , Dm = Y1 I^Trf/^ - ^Tr[/^| , (112) 
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as proposed in 0. The choice of metric does not matter much as long as it is gauge 
invariant; the two choices above give almost exactly the same Lyapunov exponent 
for a run, and cannot be used as independent measurements. To eliminate transients 
1 throw out all data until the metric distance has grown by a factor of 500; 1 then 
follow the divergence of the two paths, stopping when the metric distance is within a 
factor of 500 of saturation. This gives a very conservative fiducial period of evolution, 
over which the metric distance grows by a factor of 10^. (This is the advantage of 
using double precision.) 1 also use lattice sizes well larger than /9l to prevent finite 
size effects. To get the Lyapunov exponent I make a linear regression fit of In De 
or In Dm against time for the data in this fiducial period. To test the technique, I 
repeated one of the measurements of Krasnitz, using the unimproved Hamiltonian 
and = 4, N = 12; my value of Xmax = -348 ± .004 is within statistical error of the 
result presented in 

The technique thus tested, 1 repeated it for several lattice spacings, using the 
improved Hamiltonian; the results are presented in Table ^ which also contains the 
results of Krasnitz ||^ for comparison. It is very clear from the table that PiXmax 
has a good fine lattice limit, both for the improved and the unimproved Hamiltonian. 
It reaches its limit much faster in the improved case. The limit for the improved 
Hamiltonian differs from that for the KS Hamiltonian by about 5%, significant at 
about 5a, confirming that while the fine lattice limit of the infrared dynamics exists 
within a given lattice regulation, its value depends on the specifics of that regulation. 

As a final test of stepsize dependence, 1 recomputed the Pl = 4 datapoint using 
At = 0.1 and 0.20 (for both the thermalization algorithm and the evolution); the 
results for PiXmax are 1.147 ± .008 and 1.200 ± .010, respectively. This agrees to 
within error with a quadratic dependence in stepsize. Back-extrapolating to zero 
stepsize changes the results in the table by less than statistical error. 

Note also that the limit of the Lyapunov exponent in both regulations falls well 



below twice the damping rate of a plasmon at rest, which is [16 



2^^^5^ = .3520/r = l,408«!l. (113) 



247r " 4 

which would correspond to a value for ^L^max of 1.408 



5 Conclusion 

An improved Hamiltonian can be constructed for the evolution of the classical Yang- 
Mills equations of motion on a lattice. Using it to compute the maximal Lyapunov 
exponent, I find a much faster approach to the continuum limit, but a value of that 
limit which differs by about 5% from the value in the Kogut-Susskind Hamiltonian 
system. For the rate of Nqs diffusion, the convergence to a fine lattice limit is 
apprarently no faster, and the answers agree to within error, but the statistical power 
of the test is limited, leaving open the possibility that, as in the case of the maximal 
Lyapunov exponent, there are lattice regulation artefacts which do not vanish in the 
a = limit. 
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Such an artefact presumably arises through interactions between the infrared 
modes, which determine the quanitites of interest, and the ultraviolet excitations, 
which move under lattice dispersion relations which differ between the two regula- 
tions, and are in neither case rotationally invariant or ultrarelativistic. This possibility 
arises because the "hard thermal loop" effects which the ultraviolet modes induce are 
linearly divergent in the classical thermal theory, so most of the contribution comes 
from the momentum region where lattice effects are significant; and because the func- 
tional form the hard thermal loop effects take is not strongly restricted by symmetries 
of the theory, so the lattice dispersion relations of the ultraviolet modes can cause 
the induced hard thermal loops to have the wrong dependence on g/go- Bodecker et. 
al. have explicitly verified that this is the case for unimproved abelian Higgs theory 



111 ; in fact they show that the induced corrections to the infrared propagator are 
not even rotationally invariant. There is every reason to believe that the same thing 
happens in lattice Yang-Mills theory, in both the improved and unimproved cases. 

What might at first seem surprising is that the overall magnitude of the hard 
thermal loop effects, which in thermal units diverges as 1/a, apparently does not affect 
the infrared dynamical properties considered here, but the functional form of the hard 
thermal loops may. This is also true for one of the few infrared dynamical properties 
which is known for the quantum theory, the damping rate of a plasmon at rest; in 
its calculation the presence of hard thermal loops is essential, and their dependence 
on (f/go determines the result, but their overall magnitude turns out to cancel [|16 



Something similar seems to be happening for the maximal Lyapunov exponent. Since 
neither Hamiltonian considered to date produces the correct hard thermal loops, we 
cannot take either limit to be correct. The fact that they differ only by a modest 
amount, and that both are correctly accounting for the thermodynamic distribution 
of initial conditions, suggests that we can take the answers determined from either 
Hamiltonian as a rough (~ 10%) estimate of the correct Lyapunov exponent and 
rate of Ncs motion. However, to improve on this accuracy in the case of the rate of 
Ncs motion it will be necessary not only to determine with greater resolution what 
the rate is in the Kogut-Susskind Hamiltonian system, but to check with the same 
resolution that it is the same in the improved system. If the answers do differ at 
the level of 5 — 10%, then to get the correct rate, some technique must be developed 
which incorporates or generates hard thermal loops of the correct form. 
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Figure 2: Rate of motion of Ncs under a chemical potential for several lattice in- 
verse temperatures (reciprocal lattice spacings). The squares are the data using the 
improved Hamiltonian, the triangles are data using the Kogut Susskind Hamiltonian, 
taken from IHl 
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